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A dynamic subgrid-scale model 
for LES of the G-equation 

By A. Bourlioux 1 , H. G. Im 2 and J. H. Ferziger 3 


1. Introduction 

Turbulent combustion is a difficult subject as it must deal with all of the issues 
found in both turbulence and combustion. (We consider only premixed flames in 
this paper, but some of the ideas can be applied to the non-premixed case.) As in 
many other fields, there are two limiting cases that are easier to deal with than the 
general case. These are the situations in which the chemical time scale is either 
much shorter or much longer than the time scale associated with the turbulence. 
We deal with the former case. In this limit, the flame is thin compared to the 
turbulence length scales and can be idealized as an infinitely thin sheet. This is 
commonly called the flamelet regime; it has been the subject of many papers and 
the basis for many models (see, e.g., Linan Sz Williams 1993). 

In the flamelet model, the local flame structure is assumed to be identical to the 
laminar flame structure; thus the flame propagates normal to itself at the laminar 
flame speed, Sl ■ This allows the use of simple approximations. For example, one 
expects the rate of consumption of fuel to be proportional to the area of the flame 
surface. This idea allowed Damkohler (1940) to propose that the wrinkled flame 
could be replaced by a smooth one which travels at the turbulent flame speed, Sr, 
defined by 

St/Sl = Ai/Ap (1) 

where Ai is the total flame surface area and A p is the area projected onto the mean 
direction of propagation. This relation can be expected to be valid when the flame 
structure is modified only slightly by the turbulence. A measure of the degree of 
modification is the Karlovitz number, A'a; Eq. (1) should hold when this parameter 
is not too large. 

More recent approaches have attempted to relate the turbulent flame speed to 
turbulence intensity, it', which, presumably, characterizes the wrinkling of the flame. 
These result in relationships that typically take the form: 

S T /S L = l + C(u'/S L ) a . (2) 

For the turbulent flow dominated by an inertial range, Pocheau (1992) derived a 
linear relation (Eq. (2) with a = 1). Earlier work (Clavin & Williams 1979) also 

1 CERCA, Universite de Montreal 

2 Center for Turbulence Research 

3 Stanford University 


138 


A. Bourlioux, H. G. Im & J. H. Ferziger 


predicted linear behavior in the limit of high turbulence intensity; other authors 
have produced theories that give different exponents (Yakhot 1988, Kerstein Sc 
Ashurst 1992). In §3, we shall use DNS to demonstrate that the linear relation is 
valid, at least in the limit appropriate to LES; we shall use the database of Trouve 
Sc Poinsot (1994) and the zero heat-release data of Im (1996). 

2. Large eddy simulation based on the G-equation 

The above ideas can be applied to modeling the small scales of the wrinkling and 
thus produce the basis for large eddy simulation. It can be shown that a surface 
that moves at speed Si normal to itself in a moving fluid can be represented as a 
level surface of the G-equation (Kerstein et ai 1988): 

w + £- ( “’ G > = s ^- < 3 > 

To perform large-eddy simulation, the G-equation is filtered to produce an equa- 
tion for G, a quantity that is smoother than G. This equation contains, of course, 
terms representing effects of the scales that have been filtered out; these are the 
subgrid scale terms that must be represented by a model (Im 1995). In a simulation 
based on the filtered G-equation, the propagating flame is considered a contour of 
G, which must propagate at a speed, S, greater than the laminar flame speed; the 
increased speed plays the role of a subgrid scale model; alternative approaches to 
modeling will be discussed later. The filtered G-equation can then be written: 

W + Wj ( “ jG) = * |VG| ' (4) 

where 

S/Si = 1 + C(u 9 /Sl), (5) 

and u 9 is the velocity fluctuation characterizing the unresolved scales and C is a 
constant that can be prescribed or calculated by the dynamic procedure. In an LES, 
u f must be modeled as well. However, as the present study is an initial investigation 
of modeling the G-equation, we shall calculate u' directly from the DNS velocity 
field. Likewise, the filtered velocity field, u t is obtained directly from the DNS field. 
For later reference, we note that when the flame stretch and the Karlovitz number 
are not negligible, a diffusion-like term that represents the effect of stretch (Matalon 
& Matkowsky 1982) on flame propagation should be included on the RHS of Eq. 
(3) or (4). 

3. A priori test of the flame area scaling law 

We now present an a priori test of a dynamic subgrid-scale model for the turbu- 
lent flame speed; it is based on the model introduced by Bourlioux et ai (1996). 
Combining Eq. (1) and Eq. (5), we obtain the following equation: 

S/SL=A L /A = l + C{u'/S L h 


( 6 ) 


A dynamic SGS model for the G-equation 


139 


where Al is the flame area computed in a DNS, A is the filtered flame area to be 
computed in an LES, and 5 is the flame speed used in the LES. The latter was 
discussed above and should be selected to guarantee the correct overall burning 
rate, i.e. SA = SlAl . Eq. (6) is a useful subgrid-scale model if one can specify 
the model parameters appropriately. In this section, we first validate the linear 
relation, Eq. (6), and determine the constant C by using the resolved flame area 
Al computed from the DNS results. 

SA Test procedure 

To check the validity of Eq. (6), we process the DNS databases in the following 
way. 

1. Identify a flame surface in the DNS G-field and compute its area Al by triangu- 
lation. 

2. Filter the G and u obtained from the DNS database at various filter sizes (2A, 
4A, . . . are used, where A is the DNS mesh size). 

3. For each filter, compute u' as the square root of the subgrid kinetic energy (i.e. 
the L 2 norm of the difference between the resolved DNS velocity field and the 
filtered field). 

4. Given the filtered G-field, identify the ‘filtered flame surface’ and compute its 
area .4. 

We then investigate the relationship between the ratio Al/A and u r as a function 
of filter size. 


S.2 Databases 

We first consider premixed flames in homogeneous decaying turbulence (Trouve 
&; Poinsot 1994). In the flame, there is a smooth transition of the reaction progress 
variable, c, from 0 (fresh mixture) to 1 (burnt gas). We begin by defining a flame 
surface. Following Trouve L Poinsot (1994), we choose the level surface with c = 0.8 
as the flame front and define G = c — 0.8. Heat release effects are included in this 
DNS. Since the viscosity depends strongly on temperature, the turbulence intensity 
varies significantly across the flame so one must be careful when computing the 
turbulence intensity u'\ the value on the unburnt side of the flame should be used. 
In practice, we obtain v! by taking a 2-D Fourier transform of the velocity field on 
cross-sections ahead of the flame, averaging over the unburnt side of the domain. 
Our tests show that the choice of averaging volume is not important as long as it 
is sufficiently far from the boundary. 

The second data set is the result of DNS of the passive G-equation in forced 
isotropic turbulence (Im 1996). There is no heat release in this simulation. Several 
simplifications are used in the test procedure: 

1. The G variable is available from the DNS and can be filtered for a prion tests. 

2. In the absence of heat release, any contour of G can be considered a flame sur- 
face. The average front area can be computed from the volume average of |VG| 
(Kerstein et al. 1988). 
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FIGURE 1. Ratio of the resolved DNS flame area to the filtered flame area as a 
function of the subgrid kinetic energy u' . DNS data by Trouve &: Poinsot (1994). 

3. The entire flow field can be used to estimate u f - one does not need to distinguish 

the ‘burnt’ and "unburnt 1 regions. 

8.3 Results of the a priori test 

Results of a priori tests applied to the database of Trouve et ai (1994) are shown 

in Fig. 1. The ratio of the DNS flame area Ai to the filtered flame area A is plotted 
vs. the subgrid kinetic intensity u 9 at various times; the times shown are normalized 
by the large eddy turnover time. The DNS data were obtained on a 128 3 * * * * * grid; each 
circle is a data point; the field was filtered to grids of 64 3 , 32 3 , 16 3 (F-16 on the 
figure), 8 3 (F-8) and 4 3 (F-4). 

Figure 1 clearly shows that, if Eq. (6) is valid, its coefficient is strongly time 
dependent. There are two reasons for this. Firstly, the flame is initially planar and 
a few eddy turnover times are required to reach an ‘equilibrium’ state. Secondly, the 
turbulence is not forced; its decay can be seen from the decrease of the turbulence 
intensity at the coarsest filter size (F-4) with time. Nevertheless, the results do seem 
to indicate the existence of a universal relationship after the flame is sufficiently 
wrinkled; the change between times t=2.4 and t=4 is small compared to the change 
from t=0.4 to t=2.4. Even at large times, a distinction must be made between the 
behavior at small scales (the 64 3 , 32 3 , and 16 3 filters) and the large scales. The 
linear fit (6) appears reasonable for the small scales but not the large ones. This is 
an argument in favor of LES; modeling may be more universal for the small scales 
than for the large scales. 

In the passive database (Im 1996), the flame front is again initially planar but, 
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after several turnover times, the flame area levels out. Plots of the inverse of the 
filtered flame area vs. the turbulent intensity are very similar to those found from 
the Trouve/Poinsot database, with linear behavior for small values and quadratic 
behavior at larger scales. 

We next test the dynamic procedure; it is based on the dynamic model for noil- 
reactive flows. The parameter is adjusted using the smallest resolved scales of 
an LES. We shall not address the question of estimating u ' but focus instead on 
estimating the subgrid flame wrinkling. Given A\ and A 2 , the flame area at filter 
sizes Ai and A 2 , and the corresponding subgrid turbulence intensities, i/j and u f 2 , 
we use Eq. (6) to obtain: 

Ai/ A\ = 1 -f Cu \ , 

Aij A 2 = 1 + Cv>2‘ 

This system can be solved to produce the resolved flame area Ai and/or the model 
constant C dynamically. Table 1 gives the results for the flame speed (characterized 
by l/Ai) obtained by applying the procedure described above to Im’s data. The 
DNS field was filtered to 32 3 (F\) and 16 3 (F 2 ) grids. The modeled turbulent speed 
is compared to the exact value obtained from the DNS. The agreement is excellent. 
There is little wrinkling on the small scales and the enhancement of the flame speed 
(the difference between turbulent and laminar speeds) is very small. The table also 
gives the error in the enhancement, which is acceptable. In the next section, we 
will describe attempts to incorporate this procedure into a dynamic LES. 


Time 

7.77 

8.08 

8.37 

8.68 

8.99 

9.31 

Sdyn 

1.0881 

1.0598 

1.0825 

1.0887 

1.1246 

1.0706 

Sexact 

1.0906 

1.0633 

1.0804 

1.0801 

1.1069 

1.0672 

Ei (turbulent speed) 

-0.7% 

-0.3% 

0.2% 

0.8% 

1.6 % 

0.33 % 

E 2 (enhanced speed) 

-8% 

-6% 

3% 

11% 

16% 

5% 


Table 1. A priori test of a dynamical model for S. 

4. LES modeling test with spectral method 

In this section, the SGS model presented in §2 is tested by applying it to flames 
in forced three-dimensional incompressible homogeneous isotropic turbulence; the 
simulations are carried out using a spectral method. Heat release is neglected in 
this test so the G'-field behaves essentially as a passive scalar. 

The calculation procedure is as follows. The flow field is fully resolved on a 64 3 
grid using a pseudo-spectral method and second-order Runge-Kutta time-stepping 
(Rogallo 1981). The Reynolds number based on Taylor microscale is about 74. 
The turbulence is forced at the lowest wavenumber to maintain the kinetic energy 
constant. At every time step, the flow field is filtered onto a 32 3 grid; the resulting 
velocity field is then used in solving the filtered G-equation (4). 
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FIGURE 2. Spectra of turbulent kinetic energy and scalar fluctuations in the 64 3 
DNS calculation with u r /Sl = 0.5. 

When the code was run with the model described above, numerical instability 
resulted. Investigation showed that the instability is due to an increase in the high 
wavenumber G-field, i.e., to cusp and strong gradient formation. It is necessary to 
do something to stabilize the calculation; one possibility is to add a second order 
diffusive term, W 2 G , to the RHS of Eq. (4). As mentioned earler, similar terms 
are used to represent the effects of flame stretch and curvature on the flame speed. 
In the present DNS, T>/v — 4 is used where v = 0.015 is the molecular diffusivity. 

Figure 2 shows the spectrum of turbulence kinetic energy and the scalar fluctu- 
ations, ( g 2 ) = ((G — G) 2 ) for u'/Sl = 0.5 obtained from the DNS. The turbulence 
was forced to allow attainment of steady state spectra in a few eddy-turnover times. 
We shall use this DNS field to construct the initial condition for the LES. 
Specifically, the following subgrid-scale models are tested: 

A. S — Sl, i.e. no subgrid-scale model is used. 

B. 5 = 1.085l, where the constant 1.08 was obtained from the a priori test. 

C. S/Sl = 1 + 0.41 1( u'/Sl )> a curve fit obtained from the a priori test similar to 
Fig. 1; it' is computed from the DNS flow field. 

D. S/Sl = 1 + C(u f /S l) with the parameter C computed dynamically by filtering 
the G-field to 16 3 resolution and assuming that the model S/Sl = 1 + C(u'/Sl) 

applies at that level. The ratio S/S can be computed as the area of the constant 
G surface at the appropriate level of filtering. 

The predictions produced by all these models are compared with results obtained 
by filtering the 64 3 DNS G-field. The Markstein diffusivity V used in all of the above 
LES was increased to twice the value used in the DNS to achieve stability. This 
can be interpreted as an extra subgrid-scale scalar transport required to represent 
the effects of the filtering. A more rigorous treatment of this term is necessary 
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FIGURE 3. The contours of G obtained with various LES models after one eddy- 
turnover time; (a) G = 1.45 which is the average value and (b) G = 1.95. A single 

slice in 3-D is shown. In each figure, DNS result: ; model A: ; model 

B: ; model D: . Flame propagates from right to left. 

in the future; for example, a dynamic computation of T> can be appended to the 
computation of C . 

Figure 3 shows two G contours (1.45, 1.95) after one eddy turnover time for the 
various LES modeling strategies. Although the four fronts approximately repro- 
duce the smoothed DNS contour, the various models give different average flame 
locations. It appears that model B overpredicts the turbulent flame speed, while 
model A underpredicts it, as expected. 

The volume- averaged front location predicted by each LES model is compared to 
the DNS result in Fig. 4. All three LES models (B,C,D) overpredict the average 
flame speed. However, gradual improvement is obtained as the level of complexity 
of the model changes from the simple a priori procedure to the more sophisticated 
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Figure 4. The volume- averaged front location according to DNS and various 
LES models; model A: ; model B: ; model C: ; model D: . 


dynamic model. It is interesting to observe that, although model A underpredicts 
the flame speed, it gives the best agreement with the DNS. 

It should be noted that, the models used here are not complete. The supple- 
mentary subgrid-scale transport in the LES simulation was imposed by ad hoc 
adjustment of the diffusion coefficient, and was chosen mainly to achieve numerical 
stability of the spectral method. Improvement might be obtained by introducing a 
Smagorinsky-tvpe model for the subgrid transport with model constant determined 
dynamically. This is currently being investigated. 

5. Dynamic LES using a high order upwind scheme 

The numerical stability issue addressed in §4 is easily understood by examination 
of Fig. 5 which gives several contours G at t = 0.25. The major cause of instability is 
the formation of cusps, which are present even at this early stage of the computation. 
Another difficulty is the squeezing together of contours, resulting in high gradients 
that are difficult to capture numerically. To address those difficulties, we repeated 
some of the experiments of §5 using a different solver for the G- equation, while 
retaining the spectral velocity field computation. The G-equation is now solved 
using a numerical strategy based on level-set technology (Osher &: Sethian 1988, 
Sussman et al. 1994; for combustion applications Zhu &: Sethian 1992, Klein 1995). 
For the advection term, we use a higher order upwind code developed by LeVeque 
(1993). The source term on the right hand side of Eq. (1) is solved with the 
procedure of Zhu & Sethian (1994). A reinitialization procedure is performed at 
every time step; the G-function is reinitialized to be the signed distance function 
with respect to the flame, using the procedure of Sussman et al. (1994). This means 
that only the G = 0 contour is considered to be a flame. Figure 6 displays contours 
obtained with this method; accuracy is maintained even when the flame becomes 
very distorted and no additional numerical viscosity needs to be added when the 




A dynamic SGS model for the G-equation 


145 



FIGURE 5. Contours of the G function computed with spectral method. Flame 
propagates from right to left. 



Figure 6. Contours of the G function computed with higher order upwind 
method. Flame propagates from right to left. 

mesh is coarsened. This method is roughly equivalent to introducing a viscosity or 
diffusivity selectively at those points at which the method of the preceding section 
had trouble, i.e. at cusps and in regions of large gradient of G . 

In Fig. 7, we compare the results of a 32 3 computation with this procedure with 
the fully resolved 64 3 case results. The turbulent flame area is plotted as a function 
of time for different resolutions. The solid line is the 32 3 LES result, the dot-dash 
curve is the flame area in the 16 3 G-field obtained by filtering the 32 3 G-field. 
Using those data and the dynamic procedure of §3.3, the wrinkled flame area is 

extrapolated to the 64 3 grid ( ). This result is compared to the wrinkled flame 

area computed directly on a 64 3 grid (o ). It is clear that the suggested procedure 
underestimates the flame wrinkling. To explain the difference, we also show the 
flame area on a 32 3 grid (□ ) and a 16 3 grid (v ) obtained by filtering the DNS data. 
The discrepancy between the LES and DNS results can be traced to two effects: 
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FIGURE 7. Turbulent flame speed as a function of time from various DNS/LES 
calculations obtained with higher order upwind scheme. DNS results: 64* , o ; 32 3 , 
□ ; 16 3 , v . LES results: 64 3 , ; 32 3 , ; 16 3 , . 

- Underestimation of the flame area on the 32 3 grid: by comparing the 32 3 LES 

( ) and 32 3 filtered DNS (□ ) flame areas, it is clear that the upwind/reinitiali- 

zation scheme smooths the wrinkled front slightly. This is to be expected from 
an upwind method - this effect is relatively small and could be controlled with a 
higher order method or a solution-adaptive integration procedure. 

- Poor extrapolation of the subgrid wrinkling: extrapolation from the 32 3 and the 
16 3 grids to the 64 3 grid magnifies the error which is relatively small on coarser 
grids; this is the major source of error. 

This error can be better understood by looking at Fig. 8. In the computations, 
the linear fit (6) was used for dynamic extrapolation — the plot in Fig. 8 indicates 
that, in the early stages of flame wrinkling, the linear fit is inappropriate; a square 
root fit would be more suitable. This is consistent with the a priori test results 
reported in §3. Longer computations are being performed to assess whether this 
effect will disappear as the flame becomes sufficiently wrinkled. 

6. Conclusions 

Large eddy simulation will be necessary if reacting flows in complex geometries 
are to be simulated. This paper is a first attempt at evaluating models of subgrid 
scale effects that could be used in those flows. The laminar flamelet regime is 
considered in this paper such that the G-equation can be used as the basis for the 
modeling. 

Since the effect of filtering is to smooth a wrinkled flame, a natural model is one 
in which the smoothed flame has a higher speed than that of the laminar flame. 
Simple models of this kind were constructed and tested using the a priori approach 
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FIGURE 8. Ratio of the resolved DNS flame area to the filtered flame area as a 
function of the subgrid kinetic energy u' at various filter levels. Results for passive 
G-field with higher order upwind method. 

and large eddy simulations. A priori tests show that a linear relationship between 
the flame speed and the subgrid scale turbulent velocity is reasonable. 

The models were then tested in two types of LFS. In the first, the passive G- 
equation is solved along with the Navier-Stokes equations using a pseudo-spectral 
method. This approach is incapable of allowing heat release. Several versions of 
the model for the G-field were used including ones with a fixed constant and others 
with the parameter computed dynamically. These computations are numerically 
unstable, a problem that can be traced to the creation of cusps and high gradient 
regions. This problem can be eliminated through the addition of a diffusive term to 
the subgrid scale model. This can be justified in the same way that the Smagorinskj 
model is justified but, in this paper, the addition of the diffusive term was done in 
an ad hoc manner. 

In the other type of LES, the G-equation is solved using a high order upwind 
method and the G-field is reinitialized at each time step. This approach essentially 
introduces diffusion where required to prevent the formation of cusps and high- 
gradient regions and requires no explicit diffusive terms. 

The results show that the models are reasonable, but it appears that the LES 
models either overestimate (with a spectral method) or underestimate (with an 
upwind method) the turbulent flame speed. The reasons for this behavior are under 
investigation. 
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